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\Q ■ The effect of diffusively correlated spatial fluctuations on the proliferation-extinction transition 

of autocatalytic agents is investigated numerically. Reactants adaptation to spatio-temporal active 
regions is shown to lead to proliferation even if the mean field rate equations predict extinction, in 
agreement with previous theoretical predictions. While in the proliferation phase the system admits 
j^JTJ' a typical time scale that dictates the exponential growth, the extinction times distribution obeys a 

' power law at the parameter region considered. 

^ . PACS numbers: 05.40.-a, 64.60.Ak, 64.60.Ht 
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I. INTRODUCTION 

O ' Catalyst induced proliferation systems are very common in a variety of fields in nature such as biology, chemistry, 

S physics, and even in finance and social sciences. Technically, there is one type of agents that catalyze a reaction 
among other species, like enzymes in living animals, catalysts in chemical reactions or traders in social networks. In 
+Jj \ many cases the catalyst is not influenced by the reaction it facilitates. The spatio temporal dynamics of the catalysts 
agents dictates the dynamics of the reactive system. The theoretical understanding of such processes is usually based 
on partial differential equations for population densities, like the rate equations or (if the system is not well mixed) 
reaction-diffusion equations. Although very useful, these equations neglect the intrinsic noise build up in any realistic 
system due to the stochastic motion of individual reactants. Intuitively, the approximation involved in the study of 
deterministic processes for continuum variables rather then the noisy wandering of discrete agents is valid at high 
densities, fast diffusion limit. Unfortunately, no general systematic methodology for consistent perturbative expansion 
around this "mean field" limit is available, and the corrections are calculated on the basis of an ad- hoc methods for 
spacial cases pj or by numerical simulations . 

In a model presented recently 0, 0] the disagreement between the deterministic rate equations and the realistic 
stochastic process is emphasized for a very simple and generic system. The model includes two species: an immortal 
catalyst A that only diffuses randomly in space, and a reactant agent B which decays at rate /i and proliferate in the 
, presence of A-s at rate XNa, where Na is the number of A agents (local density of A) at the reactant spatial location. 
■ Both A and B undergo diffusion with rates Da and Db respectively. Schematically, the local reactions considered 
0^ \ are: 

<N 

qp ■ b — > (l) 

B + A — ► B + B + A. 
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The continuum description for this process is given by the mean- field rate equations for the densities of A and B, 
a(x,t) and b(x,t), respectively. The corresponding reaction diffusion equations, 

\: ' = D A V 2 a(x,t) 

■ m 

= D B V*b(x, t) - - Xa(x, t)] b(x, t), (2) 

O ■ 

admit a simple solution: since A only diffuses, catalyst density becomes a constant tia (where ua is the average A 
density) after long time, and the second rate equation becomes linear: 

^2H = D B V 2 b(x,t)-mb(x,t), (3) 

where m = fi — Xua is the decay/growth rate of the system, depending on its sign. This mean field theory predicts 
a phase transition at m — 0. For positive m the reactant concentration decays exponentially, while negative m 
yields exponential growth (proliferation). In realistic systems, of course, one expects some saturation mechanism that 
prevents explosion, perhaps in the form of B agents competition for resources. If this work this process is neglected 
and the system admits two fixed points, 6 = and b — oo. The significance and limitations of this approximation are 
discussed in the last section. 

In previous work 0, the effects of stochasticity due to the discrete character of the reactants was considered for 
this "AB model". The main observation is that diffusion never kills spatial fluctuations for discrete A agents: instead, 
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one gets at any time Poissonian fluctuations of the A density, entailing that at any time there will be a finite chance 
for a site to be "active" (with local m(x) = fi — XNa(x) smaller than zero) even if the average m is positive. At each 
of these " oases" the B population grows exponentially and a colony of B agents is developed. This leads to a shift of 
the transition from m = to finite, positive m. 

Let us present the main results of |3|, |j] . The technical tool used is a general approach for the consideration of the 
effect of stochastic noise in reactive systems on spatial domains, suggested and applied by various authors 0,0,0)- 
This theoretical framework is based on the exact Master equation of the full stochastic process. Translating the 
Master equation to the form of Schrodinger like equation for second quantized fields one may find the effective action 
and integrate fast modes to get the same effective action with reaction parameters that depend on time and length 
scale (renormalization group flows). Here, the RG transformation involves the rescaling: 



x 



t -> s z t, 



// — x </. (4) 



b -> s-%, 
A -» A/s, 

where s is the scaling factor (s > 1), d is the dimension and A is the upper momentum cutoff. This rescaling yields 
the flow equations for the model parameters: 

dX ,„ AA A 2 A d - 2 
(2 - d)X + 



dins v ' 2irDl + ^j 
dm _ n A X 2 A d ~ 2 

dins ~ m 2nDl + ^ ( ' 



driA 
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where the average diffusion rate is D = (D a + Db)/2. These flow equations were obtained by one loop perturbative 
expansion of the effective action around its mean field limit, where X/D and m/D are the small parameters. One can 
see that there is a negative correction to m and this implies that the effective "mass" m may flow to the proliferation 
phase m < as the system size increased, even if the bare mass is larger than zero. A sketch of the flow lines according 
to Eqs. J5J is presented in Figure Q. 

The aim of this work is twofold: first, to compare the theoretical predictions about the location of the extinction 
transition with numerical simulations; second, to study the nature of that phase transition. Clearly, the existence 
of proliferation phase with diverging number of reactants compels a simulation that involves with large number of 
agents, and is possible only when D is comparable with the reaction parameters. Accordingly, the RG results should 
be taken only as a qualitative guides, as the "small parameter" of the expansion is 0(1) in the numerical experiments. 
In Figure (JIJ these theoretical predictions for the extinction-proliferation phase transition, along with few RG flow 
lines are shown for d = 1, Ua = 1, D = 2 and A = 2ir /oo where ao — 1 is the lattice constant. 

The intuition gathered from the RG approach is that there are two parameter regions. In the proliferation phase 
the effective coupling constant A grows with the basic length scale of the system and for infinite sample this process 
leads to negative m. In the extinction phase, on the other hand, the growth of m dictates the system and the flow is 
toward extinction. 

In order to implement this intuition to finite samples, one should bear in mind that positive mass corresponds 
to extinction. Thus, the crossover to proliferation depends on the ability of the system to support more and more 
rescaling transformations: even if the bare mass is in the proliferation region of Figure ^ system of finite size may be 
in the extinction phase, depending on the bare parameters A, m and ua- 

In the next section the numerical technique is presented, the location of the transition is identified for finite samples 
and the results are compared with the theoretical predictions. The third section deals with the statistical properties 
of the results in different regimes where the essential difference between the phases is emphasized. In the last section 
we discuss the results and compare our study with known facts about extinction transition for bounded growth like 
the contact process in homogenous and heterogenous environments. 
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FIG. 1: RG results for AB model phase transition, Da = Db = 2, d = 1, ua = 1 and A = 2n/ao where ao = 1 is the lattice 
constant. Above the dashed line is the extinction phase for infinite system while below is the proliferation region where the 
death rate m flows to negative values. One can see that catalyst fluctuation induce proliferation in area where the mean-field 
predicts extinction (m > 0). Finite system size implies truncation of the RG flow lines, so the final effective m value may be 
positive even if the bare mass is in the proliferation region. 

II. NUMERICAL STUDY OF THE EXTINCTION TRANSITION 
A. Simulation technique 

In order to simulate the AB model, a one dimensional array of L lattice points with periodic boundary conditions 
is used. Each site i is associated with an integer (N\) number of catalysts and another integer (N B ) stands for the 
number of reactants. The hopping rates of the catalyst and the reactants are Da and Db respectively. The death 
rate of a reactant is fi, while the proliferation rate at each site is given by XN\. The initial densities of the catalyst 
and the reactant are ha and ns correspondingly For each system realization the catalysts are initially randomly 
distributed. The mean-field effective mortal rate (also defined as the reactant "mass" in the corresponding effective 
field theory) is therefore m = (J, — Xtia ■ 

The main technical obstacle for a single agent based simulation of that system is the large number of reactants in 
the proliferation phase. It turns out that an event driven algorithm based on individual agents leads to diverging 
simulation times way before any conclusive statement about the actual phase of the system is extracted. In particular, 
even a case of very large reactant " colony" may be recognized as a transient that ends up with complete extinction of 
the colony. Accordingly, a site based simulation has been used. For heavy populated sites (N g > 10 6 ) the mean-field 
equations © are solved with gaussian noise proportional to the square root of the reactant density while if N B < 10 6 
a "global" Monte Carlo method is used. For example, the number of particles that decay ("died") at certain site per 
small unit time dt is determined by generating a random number taken from a binomial distribution with probability 
jidt and N B Bernoulli trials (instead of checking N B times whether a reactant died during this time step). This 
combination of two methods allow us to reach a huge number of reactants, up to Nb = 10 300 , without neglecting the 
effect of discretization at low densities. 

The simulation is initiated with reactant density ub — 100 and catalyst density ua — 1. All the results presented 
here correspond to reactant and catalyst diffusion rates Db = Da = 2. The " mass" is therefore to = fi— Xua = /x — A. 
The simulation of an individual process ends when the system reaches either extinction or proliferation. 

An extinction of the process is declared if there are no reactants in the system. The definition of proliferation phase, 
on the other hand, is more subtle as there is no upper bound for the reactant density. We thus set an upper threshold 
of Nb = 10 50 ; above this number the simulation stops and the phase of the system is identified as proliferation phase. 
This criterion was justified by increasing the upper bound to be so large that transient effects are avoided. In fact, 
setting the threshold to be Nb = 10 300 no significant differences were identified. 
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FIG. 2: Survival probability (average over 1,000 realizations) vs m, X = 1, ua = 1, nj = 100. According to the reaction 
diffusion equations the extinction phase is above m = 0. Here one finds a region of almost sure proliferation and a region of 
almost sure extinction, separated by an uncertain region. The proliferation region grows with the system size L and saturates 
at a critical value, m c , that correspond to the extinction transition for an infinite system. 



B. The transition: theory vs. simulations 



As the system is stochastic in nature, the flow to the extinction or proliferation phase is non-deterministic, and 
for the same set of parameters one may find microscopic realizations that leads to different states. In Figure J2Jl, for 
example, the probability Poo for the system to reach the proliferation phase is plotted against m for L = 128 and 
A = 1, where each point reflects the fraction of systems (out of 1000 runs) that reach proliferation. According to 
the mean field theory © above m — the system should flow to the extinction phase. In contrary, here one finds a 
region of almost sure proliferation and a region of almost sure extinction, separated by an uncertain region. In what 
follows (for the sake of comparison between the numerics and the analytic prediction for the location of the phase 
transition) the definition of the extinction transition is taken at m for which Poo = 0.5, in order avoid the stronger 
effects of noise as P^ approaches zero or one. 

Using that analysis it is possible to continue our quest of finding the critical mass m c for the (infinite size system) 
extinction-proliferation phase transition. As a first step, the system size L is multiplied by 2 while the reaction rates 
and agents densities are kept constant. Next the survival probability Poo is plotted vs m for each system size (in 
order to change m only fj, is modified while A and ua are kept constants). As mentioned, the phase transition for a 
finite size system is defined where Poo = 0.5. The critical m of the finite size system is plotted versus InL in Figure 
(PJa). One can see that for small system sizes, the critical "mass" increases as the system size L increases, since the 
probability for larger catalyst fluctuation grows with size. There is, however, an intrinsic limit (independent of the 
system size) on these spatial fluctuations, and according to the RG analysis one expects that for large enough bare 
mass there will be no further increase in the critical mass as a function of the system size [see Figure j5|b)]. Indeed, 
the simulation indicates that for larger system sizes the critical mass begins to saturate while it approaches the "real" 
phase transition. Since the "real" phase transition at the thermodynamic limit is beyond our computational abilities, 
the critical mass m c for certain A is evaluated by fitting the results to the function 

m c 

m(x) = m c - 6 
1 + (x/x ) p 

where x = ln(L). As indicated in Figure 10b), this function yields an almost perfect fit to the RG based theory, and 
seems also to fit the numerical experiment quite well [Figure (0 a)]. 

After the critical mass for a specific A is determined, using the above procedure, one proceeds to plot the infinite 
system size phase transition line m c vs A, Figure Above the transition line is the extinction phase for the infinite 
size system, while below is the proliferation phase. The theoretical phase transition line, based on the RG analysis, 
is plotted together with the line obtained from the numerical experiment on a logarithmic scale. Both theoretical 
and numerical results suggests that m c ~ A Q , where the values of a w 1.4 (for the RG phase transition curve) and 
a « 1.45 (for the numerics) are close and the curves are almost parallel. This, again, implies a qualitative similarity 
between the numerical simulations and the picture emerges from the theory. 
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FIG. 3: m Vs ln(L): (a) Simulation results (squares) for A = 1, ua = 1, ub = 100. One can see that the system robustness is 
increasing (proliferation for larger death rates) as the system size grows. The critical "mass" m c (the mass for which there is 
no proliferation even in the limit L — > oo) is determined using the best fit to the logistic function JSJ (full line) . (b) RG results 
(squares) for A = 1, ua = 1 and A = 2n/ao where ao = 1 is the lattice constant. There is a region where the phase transition 
depends on the system size, but above a critical mass m c there is no growth independent of L, in agreement with Figure Q. 
The full line is the best fit to ©. The qualitative resemblance between the RG predictions to the simulation results is clear. 
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FIG. 4: Phase transition lines in the m - A plane. The full line is the theoretical predication based on the RG perturbative 
analysis while the squares are the results of the numerical experiment (see text). 

III. SAMPLE TO SAMPLE FLUCTUATIONS AND THE PHASE TRANSITION 

Let us take a deeper look at the dynamics of the system, trying to describe the different regions of Figure J5J. 
For a point in the distinguished proliferation region (to < 0.6, Poo ~ 1), all the samples proliferate. The growth of 
the B population is exponential, and both the average over the logarithm of the population and the logarithm of the 
average grow linearly in time [Figure JSJ)]. (if most of the realization decay and only few show exponential growth, 
(Nb) is dominated by rare events and grow, while (ln(iVB)) decays to zero since the logarithm kills the influence of 
rare events on the average.) If we further plot a histogram of the "end times" for various realization (i.e., the time 
needed to the system to reach the proliferation upper bound) one can see that there is a typical end time and the 
histogram corresponds to Gaussian distribution around its value [Figure @]. The lines of Figure JSJ do not coincide 



6 



700- 




time 



FIG. 5: \si{Nb) and Qxi(Nb)) vs time (average over 1,000 realizations), A = 1, m = 0.6, ha = 1, ns = 100. The typical case is 
shown to be similar to the average over all realizations. The lines do not collapse since (ln(A r s)) is dominated by the peak of 
the gaussian of Figure ©, while \xi(Nb) is dominated the tail of this histogram. 
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FIG. 6: End times histogram (1,000 realizations), A = 1, m = 0.6, tia = 1, "-s = 100. One can see that there is a typical time 
for the system to reach the proliferation upper bound. 

since (\ii(Nb)) is dominated by the peak of the gaussian, while Iu(Nb) is dominated by its tail. 

Deep in the extinction region (e.g., m = 2, = 0) the situation is different. In this region all the samples end up 
in extinction, but the corresponding time scale (the effective decay rate, or the effective mass) shows strong sample 
to sample fluctuations. In Figure @ a histogram of death times (the inverse mass) is shown. If one assumes that, at 
long times, the reactant concentration decays like exp(— t/r), this figure implies that P(t) ~ t~ 6 where for the results 
of figure (0) 9 « 5. The strong decay of P(t) manifest itself in the similarity between the typical and the average 
cases, Figure (JSJ. 

Close to the edge of proliferation (e.g., m = 1.1, P^ « 1) the power law for the decay is much smaller. In Figure 
© the histogram of death times indicates 9 ~ 2. Accordingly, there is a lot of difference between the average (ln(TVs) 
, Figure IjlOl a)) and the typical ((\ii(Nb)), Figure l|10lb)) time evolution of the total B population. While the typical 
case decays smoothly to zero, on average the system grows, saturates, and decays after much larger time. The samples 
for which r is large are samples where the reactant population grows rapidly on time scales smaller than r and then 
decays. Accordingly, these samples dominate the average population while the typical population is determined by 
the large number of samples that admits small r Q . 

In the middle region, say, 0.6 < m < 1.1, the survival probability is between zero and one, and accordingly two 
peaks of death and proliferation times are seen in Figure l|ll|) , corresponding to the proliferating and to the extinction 
events. Again there is a Gaussian like distribution around a typical time scale associated with the proliferation 
samples, and power law distribution of time scales in the extinction samples. 

Our numerical experiment suggest, thus, that the extinction transition of the saturation free AB model takes place 
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FIG. 7: Death times histogram (log-log scale, 200,000 realizations) for A = 1, m = 2, ha = 1, ns = 100. The death times 
histogram seems to maintain its power-law tail even above the "real" extinction phase transition. However, it decay much 
faster with exponent 8 ~ 5. The rare events in that case are "too rare" and the average is dominated by the typical events, as 
emphasized in Figure |8] 
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FIG. 8: ln(iVs) and (rn(iVs)) vs time (average over 1,000 realizations), A = 1, m = 2, tia = 1, TIb = 100. One can see that the 
average case is similar to the average over the realizations. There are no significant fluctuation (within our finite realization 
number) since finding exponentially rare events is almost impossible 

in the following way: deep inside the extinction phase there is a power law distribution of decay times with 6 > 2, 
and 8 becomes smaller and smaller while decreasing m. Close to the proliferation threshold 9^2, and from now on 
there is finite fraction of the samples that proliferate to infinity. The sample to sample fluctuations decrease as the 
transition is approached, and in the "living" phase Gaussian statistics is observed. 

IV. CONCLUSIONS AND DISCUSSION 

This paper presents a numerical experiment that confirms, qualitatively, some of the theoretical predictions about 
the extinction transition in system of autocatalytic agents with diffusively correlated fluctuations. The extensive 
simulations used render a phase diagram similar to the predictions of the exact Master equation based renormalization 
group technique, even in the parameter region where the "small parameters" of the perturbative technique are 0(1). 
As expected, resources fluctuations increase the ability of species to survive, but in order to extract the potential of 
these fluctuation via the adaptation mechanism that associate reactant colonies with catalysts islands, the size of the 
system must be large enough. This finite size effect has been demonstrated, and the system is taken to the large scale 
limit in order to approach the "real" phase transition. 
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FIG. 9: Log-log scale death times histogram (200,000 realizations), A = 1, m = 1.1, n.4 = 1, tib = 100. The probability of 
realization to decay (Nb = 0) between t to t + At is shown. The data correspond to a power law distribution of decay times 
with P(r) ~ r~ e with 6 ~ 2. This implies that, to the accuracy of our numerical experiment, the average time associated with 
the decay of the system in the extinction phase diverges. 




FIG. 10: Typical (a) versus average (b) decay patterns. In the upper panel {Iii(Nb)) is plotted against time and shows a 
smooth decay, while in the lower panel Iii{Nb) shows completely different pattern with different time scales. The average is 
taken from 200,000 different realizations, for A = 1, m = 1.1, tia = 1, ns = 100. Most of the samples decay to extinction very 
fast, but the average (b) is dominated by rare events, due to strong sample to sample fluctuations. 
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FIG. 11: Histogram of death times (a) and of proliferation times (b) in the intermediate parameter region. One sees that 
the death times distribution admits a long tail corresponds to power- low [see the inset of panel (a)], while proliferation times 
distribution is sharply peaked around its typical value. 
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Another result, not predicted by the theory, is the appearance of a scale free fluctuations in the extinction phase, 
in contrary with the Gaussian distribution associated with the characteristic time in the growth phase. 

Let us compare, now, these findings with some known results about extinction transitions. A common conjecture, 
first suggested by Grassberger and Janssen is that, in the absence of conservation laws, all the extinction 

transitions belong to the directed percolation (DP) equivalence class ^lj- This proposal is based on the following 
intuition: close to the transition point the average density of the reactants is zero. In a single time step reactant may 
decay, survive or yield an offspring. These processes are the basic rules of DP, hence all transitions are at the same 
equivalence class. 

In the presence of spatial heterogeneity, on the other hand, the situation is more complicated as quenched ran- 
domness appears to be a relevant operator. It was shown by Janssen 0] that the corresponding field theory admits 
only runaway solutions and no stable fixed point is reached in the physical domain, so the transition in that case is 
not in the DP equivalence class. This analysis supports the numerical work of Moreira and Dickman who used 
Monte Carlo simulations to analyze the 2d contact process with various fractions of randomly distributed inactive 
sites. In the extinction region these authors identify the Griffith phase, with nonunivcrsal power laws. In that regime 
the local dynamics, associated with rare events, dictates the system's evolution. While the chance to find a "good" 
region (with no inactive sites) of linear size R decays exponentially with _R, its typical lifetime grows exponentially 
with size, yielding a non universal power law below the transition. 

One may a-priori suggest that, for the AB model considered here, the diffusive dynamics of the catalysts simply 
average out their effect, thus the resulting transition should be in the directed percolation equivalence class, like the 
transition on homogenous substrate. This assumption is wrong, since a diffusive heterogeneity is a relevant operator 
at the transition. One may either apply the Harris criterion or consider the naive scaling for the action associated 
with the Reggeon field theory with diffusive disorder: 

S = { altdx + Dj^*)(Vtp) + - + <f ^ + £> (V0*)(V<£) ~ + (7) 

J at at 

to see that it supports only runaway solutions, and that flows to zero upon rescaling transformation. 

Since the system seems to flows into the quenched disorder phase one may suggest that, for an unbounded growth, 
the colonies on the rare oases grow without limit, so there is no phase transition in the thermodynamic limit (unless 
the catalyst distribution and dynamical rules are such that they never allows for m < at any finite region). However, 
Kesten and Sidorovicious |l4j | proofed that, even without saturation term, there is a parameter regime that corresponds 
to the extinction phase, where all processes decay with probability one for an infinite sample. The stochastic wandering 
of the catalysts, thus, seems to control the system and to yield a possibility of extinction even in the infinite capacity 
limit. 

The emerging picture, supported by the numerical simulations presented here, is that the phase transition exists, 
with a wide (power law) distribution of extinction times (Griffith phase) below that transition. The diffusive wandering 
of the catalysts is not strong enough in order to replace the fluctuating concentration by its average and the adaptation 
of the reactant colonies to the Poissonian noise is of crucial importance. On the other hand, even if the carrying 
capacity is unbounded this wandering is strong enough to yield a finite range of parameters where the process ends 
up at extinction with probability one. 
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